A coupled map lattice model for spontaneous pore formation in anodic oxidation 
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We construct a coupled map lattice model for the spontaneous pore formation in anodic oxidation 
in two dimensions and perform numerical simulations, after we explain steady flat solutions, their 
linear stability and a single pore solution for a model of Parkhutik and Shershulsky. 
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CN 1 I- INTRODUCTION 

H 

Nanoscale pores are spontaneously created in alumina film, when anodizing aluminum is oxidized to alumina (AI2O3) 
in various acid solutions A schematic figure of the pore structure is drawn in Fig. 1(a). The anodic aluminum 
is used as aluminum pots and aluminum window sash in daily life. The alumina film works as a protection film 
against further corrosion. The configuration of pores is usually irregular, but regular arrays of pores were found in 

1 1 controlled experiments @, 0| . Similar porous structures were reported in anodic titanium Q . The mechanism of the 

pore formation in anodic oxidation has been long studied, but is not still completely understood [f|. Recently, the 
porous anodic films were studied from a view point of the instability of growing interfaces [ta-flOl. Unstable growing 
interfaces have been intensively studied in the problems of dendrites and viscous fingering The feature in the 

problem of the anodic oxidation is that there are two interfaces between Al and AI2O3 and AI2O3 and acid solution 
(electrolyte), and their interaction is important. The linear stability analysis for flat interfaces was performed by 
1 \ Thamida and Chang Q, but the model does not provide a physically justified short-wave cutoff. An elastic effect by 
*^ ■ the volume change at the oxidation can play an important role for the formation of the regular hexagonal arrays of 
pores [3|. Singh proposed a more plausible model including the elastic effect, performed the linear stability analysis, 
and derived a weakly nonlinear equation similar to the Kuramoto-Sivashinsky equation 

However, the well-developed pores were not well treated because of the difficulty of the moving boundary conditions. 
The Stefan problem with the moving boundary conditions has been studied typically in the problem of dendritic crystal 
growth. There are several models to treat well-developed complicate interface patterns. Witten and Sander proposed 
a DLA (diffusion- limited aggregation) model for fractal growth patterns [l2j]. The phase field model was proposed 
to perform numerical simulations for dendritic patterns [13]. We proposed a coupled map lattice model as a simpler 
\ simulation model of dendrites and extended it the oscillatory electrodeposition and the pore formation in activated 
carbon In this paper, we study a simple model for the spontaneous pore formation in anodic oxidation which 

was first proposed by Parkhutik and Shershulsky, and perform numerical simulations using a coupled map lattice 
model. The model docs not include the elastic effect and cannot reproduce a regular array of pores. However, the 
model can simulate a strongly nonlinear evolution toward the well-developed pores and competitive dynamics among 
the neighboring developed pores. 
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II. SIMPLE MODEL FOR ANODIC OXIDATION AND ITS LINEAR STABILITY ANALYSIS 

The anodic oxidation of aluminum creates alumina film on the surface of aluminum. Flat alumina films are stable 
in alkali solutions but they are unstable to form many pores in acid solutions such as sulfuric acid solutions. It is 
known that the thickness d of the alumina film is proportional to the voltage V of the anode and the average interval 
between neighboring pores increases with V. The detailed processes of chemical reactions are not well known but the 
chemical reaction: 



2A1 + 30H~ -> AI2O3 + 3H+ + 6e _ (1) 

occurs at the interface between Al and AI2O3, and the alumina dissolves chemically at the interface between AI2O3 
and the electrolyte. As a result of the chemical reactions, the interface patterns as shown in Fig. 1(a) move downwards, 
i.,e., in the — z-direction. The chemical reaction rates depend on the voltage at the interface. The electric current 
is induced by the chemical reactions. We study a simple dynamical model for the metal-oxide and oxide-electrolyte 
interfaces. The positions of the two interfaces are respectively expressed as Ci(x,t) and CaC*, t), whose schematic 
figure is drawn in Fig. 1(b). The electric potential is denoted as cj>. The voltage is assumed to be = U at the first 
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FIG. 1: (a) Schematic figure of porous alumina, (b) Two interfaces of aluminum-alumina and alumina-electrolyte. 



interface £i- The electric potential satisfies 



V 2 6 = 0. 



The growing velocity of the first interface is written as 

dt 



-aE„ 



(2) 



(3) 



where E e is the electric field strength, a > is a Faradaic coefficient multiplied by the conductivity of the oxide, and 
d n implies the normal differentiation. Parkhutik and Shershulsky proposed a model for the velocity of the second 
interface at (2 as 



d( 2 

— = -a exp(k d E e ) + /3 exp(k o E e ), 



(4) 



where ao, (3q, kd and k a are some coefficients for dissolution and oxidation. Here, we propose a simpler model, assuming 
that the electric field is weak near 



at 



b 2 E 2 e + Dk, 



(5) 



where 6o,&iand b 2 are expansion coefficients of the interface velocity with E e , K = (d 2 /dx 2 + d 2 /dy 2 )^2 denotes the 
curvature of the interface, and D is a coefficient proportional to the surface tension. The interface tends to become 
flat by the effect of the surface tension. The electric potential is assumed to <f> = at z — £2- 

If the flat interfaces are stable, the electric potential satisfies <j> = V(C2 — z )/((2 — Ci), and d n (j> — —d<j>/dz = V/d, 
where d = (2 — Ci is the interval between the two interfaces. The velocity d(i/dt = —vo = —aV/d and the velocity 
d&fdt = — vq — ~b — bi(V/d) + b2(V/d) 2 . The steady electric field E = V/d is determined from dQ\/dt = d^/dt as 



V _ bi-a + ^{b x - a) 2 +Abgb2 
d ~ 2b 2 



(6) 



This relation shows that the thickness d is proportional to the voltage V. If b 2 = 0, V/d — bo/ (a — 61), and 60 > 
and a > bi are necessary for stable growth. If bo — 0, V/d — (bi — a)/&2, and 62 > and b\ > a are necessary for 
stable growth. 

The stability of the flat interfaces can be investigated, assuming the perturbation of the form: 



5cj) = V {d~ z-v t)/d + 5(f) 1 q cos(qx)e- qz+x o t + 5<f) 2 cos(qx)e qz+x -' t 1 
Ci = -v t + a q cos(qx)e Xqt , ( 2 — ~v t + d + b q cos(qx)e x,lt . 

From the boundary conditions at the two interfaces, 

-{V/d)a q + 5<j>l + &<% = 0, 
~(V/d)b q + 8<t>\e-« d + S0 2 e" d = 0, 



(7) 



\ q a q + aq(8<j) q - 84%) = 0, 

(X q + Dq 2 )b q + {h + 2b 2 (V/d)}(e 



(8) 
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FIG. 2: Eigenvalue X q at bi = 0.035 and 0.02 for a = 0.015, b = 0, , b 2 = 0.04 and £> = 0.003. 




FIG. 3: (a) A single pore solution at A = 0.5, xo = 13.5, V = 2 and &i/a = 0.5. (b) An array of pore structures, (c) Thickness 
d as a function of V for A = 0.5, a;o = 13.5 and bi/a = 0.5 

The eigenvalue X q satisfies 

_ e -?d )A 2 + [{jD(? 2 + aq( y/ d ) _ b iq (V/d) - 2b 2 q{V/d) 2 }e qd - {Dq 2 - aq(V/d) + bxq(V/d) + 2b 2 q(V / ' df}e- qd ]\ q 
+aq(V/d)[{Dq 2 - hq(V/d) - 2b 2 q{V / df}e qd + {Dq 2 + hq(V/d) + 2b 2 (V/d) 2 } e - qd } = 0. (9) 

Figure 2(a) shows \ q at b x = 0.02 and 0.035 for a = 0.015, 6 = 0,b 2 = 0.04 and D = 0.003. The flat interfaces is 
stable at b\ — 0.035, but they are unstable at b\ = 0.02. This instability is analogous to the Mullins-Sekerka instability 
in crystal growth, however, the interaction between the two interfaces is important in this problem. 

III. A SINGLE PORE SOLUTION 

When the flat interfaces become unstable, many pores are spontaneously created. A single pore solution can be 
explicitly constructed for & = &2 = in two dimensions as a generalization of the Saffman- Taylor solution in a 
channel for the viscous fingering 

|l7|, 01 • Single pore solutions for more complicated conditions such as Eq. (4) were 
numerically solved by Thamida and Chang [7|. The potential satisfies (d xx + d zz )<fi = 0, and the boundary conditions 
are ad n (j) — U cos#i at the first interface and b\d n (j) — Ucos9 2 , where U is the uniform propagation velocity in 
the — z direction, and 8i, 2 are the angles between the normal directions of the two interfaces and the z direction. 
For a spatially-periodic array of pores with wavelength 2x$, 2x$ corresponds to the channel width in the Saffman- 
Taylor solution. Another variable ip is introduced, where (/) + itf> is a holomorphic function in the complex plane 
— z + ix. From the boundary conditions ad n 4> = U cos#i, bid n (f> — U cos9 2 , dip/dx = U/a at the first interface, and 
dtfj/dx = U/bi at the second interface, which yields ip = Ux/a at the first interface, and i/j = Ux/bi at the second 
interface. Furthermore, <fi = V at the first interface, <fi = at the second interface, and ip = ±Vb is assumed at the 
boundaries x — ±xq. The width of the pore is assumed to be 2\xq. The schematic figure is shown in Fig. 3(a). The 
inverse function x(4>, ip) can be expanded with the Fourier series as 

— = — +y"{A n sm(nniP/V )e~ n ^ Vo + B n sm(nniP/V Q )e n ^/ Va }. (10) 
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FIG. 4: (a) Pore structure at a — 0.042, &o = 62 = 9 = 0, bi = 0.01 and V = 5 starting from d(0) = 16. The system size is 
35 x 400. The shaded region represents oxide, (b) Time evolution of the thickness d(t). (c) Pore structure at a — 0.042, 60 = 
b 2 = g = 0, 61 = 0.01 and V = 5 starting from d(0) = 10 



At the first interface, i/j — Ux/a and <p — V, which yields 



x U V 



Y,{A n e~ n7rV/Vo + B n e™ v / y °}sm(mnP/Vo). (11) 



At the intersections of the first interface and the boundaries x — xq, tp takes Vq owing to the continuity, and therefore 
a/U = xo/Vq, and 

B n = -e 2 ™ v ' v °A n . (12) 

At x = Alo and z = 00, ip approaches Vq. It leads to UXxo/bi — Vq and A = bi/a. At the second interface, ip = Ux/b\ 
and (f> = 0, which yields 

(A ~y = V(A„ + B„) sin(n^/K))- (13) 

From the Fourier expansion (13) and the relation (12), 

2(-l)"(l-A) 

" (n7r){l-exp(-2n7rV/Vb)}' 1 J 



Using the Cauchy-Riemann relation: dx/dtp — —dz/d(f>, the first interface satisfying <j> — V is expressed as 

z V 
xo Vo 



— = -^r + Y{2A n )e- n * v/Vo cos(mrx/x ) + C, (15) 



where an integral constant G is determined as C = V/Vq + J2(— 2A n )e nnV / v o by assuming z = at x = 0. The 
second interface satisfying <j) = is expressed as 

— = V {(A„ - B n ) cos{n7rx/(Aa;o)} + C. (16) 

Figure 3(a) is an example of pore solution for A = 0.5, x a = 13.5, V = 2 and &i/a = 0.5. Figure 3(b) shows an 
array of pore structures obtained by repeating the single pore solution by wavelength 2xq. The mirror-symmetric 
spatially-periodic solution is also a special solution satisfying the boundary conditions. The pore solution by Thamida 
nd Chang seems to have a sharp cusp at the boundaries [7], however, our solution has no such singularity. Figure 
3(c) shows a relation of the thickness d of the two interfaces at x = as a function of V. A monotonic increase of d is 
observed. When the thickness d is large, the first interface becomes flat and the second interface take a form similar 
to the Saffman- Taylor solution. 



IV. NUMERICAL SIMULATIONS WITH A COUPLED MAP LATTICE MODEL 



We propose a coupled map lattice model on a square lattice corresponding to the model expressed by Eqs. (2), (3) 
and (5) to investigate the nonlinear time evolution of the anodic oxidation. The time and the space are discretized in 
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FIG. 5: (a) (b) Pore structures at (a) V = 3 and (b) 7 for a = 0.02, b Q = g = 0, bi = 0.01 and b 2 = 0.005. 5(c) Relation of the 
thickness d and V. The dashed line is a linear fitting line. 
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FIG. 6: (a) Flat pattern at b = 0, a = 0.015, bi = 0.035, b 2 = 0.04, g = 0.0003 and V = 12. (b) Pore structure at a = 0.02, b = 
0, 6i = 0.01, b 2 = 0.005, g = and V = 2. (c) Pore structure at a = 0.02, 6 = 0, 6i = 0.01, b 2 = 0.005, g = 0.005 and V = 2. 
(d) Pore structure at a = 0.045, 6 = 0.0002, 6i = 0.01, 6 2 = 0, 3 = and V = 3. 



the coupled map lattice model. The discrete Laplace equation is written as 

cf) n {i + l,j) + <f) n {i-l,j) + Mhj + 1) + ~ 1) - #„(m) = 0. (17) 

We have solved the discrete Laplace equation by iteration of a hypothetical discrete diffusion equation. 

C +1 (^j^C(^jO + ^{C(* + L^ + C(^L^ + C(^j + i) + C(^j-i)-4C(^j)}, (18) 

with = 0.2, and the iteration number of m is set to be 400 in most numerical simulations. The potential <f> is set 
to be <j) = V in the metal region and <f> — in the electrolyte region. We introduce two kinds of order parameters p 
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FIG. 7: (a)Pore structure at a = 0.02, a = 0.01, k d = 0, /3 = 0.01, k a = -1 and V = 2. (b) Pore structure at a = 0.045, a = 
0.05, k d = 0.2, /3 = 0.0498, fc D = and V = 3. (c) Pore structure at a = 0.045, q = 0.01, fc d = 1.2, /3 = 0.005, k a = 1.1 and 
V = 3 Pore structure at a = 0.02, a = 0.01, fc d = O I( O = 0.01, fc D = -1 and V = 2. 
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and q. The order parameters are p — and q = 1 in the metal region, p = 1 and q = 1 in the oxide region, and p = 1 
and q = in the electrolyte region. The order parameter changes in time only at the interface sites. In the interface 
sites between metal and oxide, p takes a value between and 1, and p increases as 

Pn+i(i,j) =Pn(iJ) +a\d n (j) n (i,j)l (19) 
where \d n <j>\ is numerically evaluated with 

M = VWi + M) - Hhj)) 2 + W - ij) - Hhj)) 2 + + 1) - W, j)) 2 + WTj - 1) - j)) 2 }/2, 

because (V0) 2 = (d n (f>) 2 + {d s (j)) 2 and the differential 9 s </> along the interface is zero. If p n+ i(i,j) goes over 1, the 
interface site changes into a metal site. In the interface sites between oxide and electrolyte, q takes a value between 
and 1, and q decreases as 

q n +i(i,j) = q n (hj) - i b Q + hE en (i,j) - b 2 E en (i,j) 2 +g{N oc (i,j) - 3)}, (20) 

where E en = \d n <j)\ is approximated at (f> n (hj) because (j) n (i,j) is in the electrolyte region, N oc (i,j) is the number 
of the electrolyte sites which locate at the nearest and the next-nearest neighbor sites of (i, j), and g is a coefficient 
corresponding to D in Eq. (5). The term N oc — 3 is used to represent the surface tension effect, because N oc is 3 for 
flat interfaces and the surface tension effect is zero for the flat interfaces. This term was used in the coupled map 
lattice model for dendritic growth including the surface tension effect [12]. If q n (i,j) goes below 0, the interface site 
changes into an electrolyte site. Figure 4(a) displays the oxide region obtained by the coupled map lattice model with 
system size 27 x 400 for a = 0.042, bo = b 2 = g = 0, b\ = 0.01 and V = 5. The solid and dashed lines are the analytic 
single pore solution obtained in the previous section. Fairly good agreement is seen. Figure 4(b) shows the time 
evolutions of the thickness d between the first and the second interfaces, when the initial value of d(0) is changed. 
When d(0) = 16, the thickness d(t) is almost constant in time, however, d(t) decreases in time for d(0) — 10. Figure 
4(a) shows a pore structure for d(0) — 16. Figure 4(c) displays the oxide region for d(0) = 10 at n = 25000. The 
steady pore solution by Eqs. (15) and (16) is not always realized as a stable solution when bo = b 2 = 0. 

If bo = and b 2 > 0, a unique stable pore solution is obtained, because too small d makes the velocity of the second 
interface slower and d becomes larger. Figures 5(a) and (b) display the oxide region obtained by the coupled map 
lattice model of 27 x 400 at V = 3 and 7 for a = 0.02, b = g = 0,h = 0.01 and b 2 = 0.005. Figure 5(c) displays 
a relation of the stationary value of the thickness d as a function of V obtained by the coupled map lattice model. 
A monotonously increasing relation is obtained. The single pore structure is similar to the Saff man- Taylor solution, 
when d is relatively large. 

In numerical simulations shown in Fig. 4 and 5, a small dent is set in the second interface around x = 15 as an 
initial condition. In Fig. 6, we show some numerical results in a larger system of 400 x 400. The initial conditions 
are flat interfaces with small random perturbation. Figure 6(a) shows stable flat interfaces for bo = 0, a = 0.015, b\ = 
0.035,^2 = 0.04,5 = 0.0003 and V — 12. The flat interfaces are stable for the parameter values as shown in §2. 
Figure 6(b) shows a pore pattern for a = 0.02, b = Q,bi = 0.01, b 2 = 0.005,3 = and V = 2 at n = 60000. Fine 
pore structures are created initially, because g = 0, and it remains on the top region. A kind of competition occurs 
among pores, and some pores grow and the others stop to grow. After all, a rather regular array of pores appears 
below z — 250. Qualitatively similar time evolution is observed in experiments. Figure 6(c) shows a pore structure 
for a = 0.02, b = 0, h = 0.01, b 2 = 0.005, g = 0.005 and V = 2 at n = 60000. By the effect of the surface tension, a 
fine structure does not appear and the average interval between neighboring pores after the initial transient is larger 
than in the case of g = 0. Furthermore, the growth velocity and the depth of the pore is smaller than in the case 
of Fig. 6(b). Figure 6(d) shows a pore structure for a = 0.045, b = 0.0002, b x = 0.01, b 2 = 0, g = and V = 3 at 
n = 17000. The top of the second interface dissolves and sinks downwards in time because bo is not zero. 

We can perform numerical simulation of a coupled map lattice model corresponding to the original model expressed 
by Eqs. (2), (3) and (4). For this model, Eq. (20) is replaced by 

?n+i(*ji) = Qn(i,j) - [a exp{k d E en (i,j)} - /3 Q exp{k E en (i,j)}}. (21) 

Figure 7(a) shows a pore structure for a = 0.02, ao — 0.01, kd — 0,/3o — 0.01, k a — 1 and V = 2 at n — 40000. 
The Taylor expansion of Eq. (21) by E e yields bo = 0, b\ = 0.01 and b 2 = 0.005 in Eq. (20), which corresponds 
to Fig. 6(b). Figure 7(b) shows a pore structure for a = 0.045, a a = 0.05, k d = 0.2, /? = 0.0498, k Q = and 
V = 3 at n = 11000. The Taylor expansion of Eq. (21) by E e yields b = 0.0002, b x = 0.01 in Eq. (20), which 
corresponds to Fig. 6(d). Similar pore structures appear in the model Eq. (21). Figure 7(c) shows a pore structure 
for a — 0.045, ao = 0.01, kd = 1.2, /3o = 0.005, k a = 1.1 and V — 3 at n = 10000. The alumina film becomes even 
thinner by dissolution, and cusp-like structures rather than pore structures appear in the second interface. 
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V. SUMMARY 



We have studied a simple model for spontaneous pore formation in anodic oxidation. It is one of typical problems 
for growing interfaces, however, this system is unique in that the interaction between two interfaces is important. The 
second interface between the oxide and the electrolyte is essentially unstable, because the interface moves toward the 
region of large electric held, and the hrst interface between the metal and the oxide is essentially stable, because the 
interface moves away from the region of large electric held. The interaction between the essentially stable hrst interface 
and the essentially unstable second interface makes a new kind of growth pattern. We have proposed a coupled map 
lattice model to understand qualitatively the time evolution of the anodic oxidation and performed several numerical 
simulations. We have reproduced some qualitative features of the anodic oxidation, such as the single pore solution, 
the dependence of the thickness d on V, and the strongly nonlinear evolution of the pore structure. 

However, our model studied in this paper is a rather simplified one. We need to improve the model to incorporate 
the better boundary conditions using the Butler- Volmcr relation, the reaction-diffusion dynamics of some chemicals 
such as H+ and OH~. Especially, we would like to take an elastic effect into the coupled model to reproduce the 
regular array of pores. We hope that a regular hexgonal array of pores might be obtained in the three-dimensional 
simulations of the generalized coupled map lattice model. 
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